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I. INTRODUCTION 


In the present Report a detailed description is given of a basic 
computational program for the evaluation of three-dimensional, super- 
sonic, inviscid, steady flows past airplanes. For the sake of simplicity, 
no imbedded shocks are considered here. The emphasis is put instead 
on how a powerful, automatic mapping technique is coupled to the fluid 
mechanical analysis in order to assure a high degree of accuracy without 
increasing the number of computational nodes beyond reasonable limits. 

Care has been taken to describe and to code each of the three 
constituents of the analysis (body geometry, mapping technique, and gas 
dynamical effects) separately, to facilitate applications to different 
geometries or substitution of the present set of unknowns and equations 
of motion by other sets. Sections V through X contain the outline of the 
code dealing with gas dynamical effects; all their statements and formulas 
are unaffected by changes in the mapping technique or mapping parameters 
or in the geometry of the airplane. All expressions related to the map- 
ping are given in Sections XI and XII. Results of computations based on 
sample geometries, and discussions will be presented in a separate 
Report. 

II. FRAMES OF REFERENCE 

The free stream flow is assumed to be uniform, with a given 
Mach number, M^. A Cartesian, orthogonal frame of reference, (x, y, t) 
is defined as having the y and t-axes in the symmetry plane of the 
vehicle, the t-axis lying along the fuselage. The unit vectors of the x, y, 
and t-axis are called I, J and K, respectively. The free stream velocity 
vector, V^, is parallel to the (y,t)-plane; the angle of attack, a, is the 
angle between and K; therefore. 
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( 1 ) 


V = V_(Jsino + KcosO.) 

CO 00 ' 


In each cross-sectional plane, a complex variable, z, is defined as 
z = X + iy (2) 

A conformal mapping (details of which will be found in Section XII) defines 
a one-to-one correspondence between the portion of interest of the right- 
hand side of the z-plane and a portion of the right-hand side of a C-pls-ne 
where, by and large, the image of the cross-section of the airplane is 
nearly circular; it is convenient, thus, to express the complex variable 
^ in the form: 


C = p e"® (3) 

The analytic function C(z) implies that p and 9 are functions of 
X and y, and vice versa. Such functions, in general, change from one 
cross-section to another; therefore, we may write: 

^ fx = x(p, 0 ,t) 


P = P(x, y, t) 
0 = 6(x, y, t) 
T = t 


y = y(p, 

t = T 


(4) 


We must take good care of denoting t by another sy7nhol, T, when con- 
sidered in connection with p and 6 since when t changes and x, y 
remain unchanged, p and 6 generally change; consequently, derivatives 
with respect to t (at constant x and y) generally differ from derivatives 
with respect to T (at constant p and 6). Let p = b(9, t) and p = c(9, T) be the 
equations of the image of the airplane body contour and of the image of the bow 
shock in the C- plane. A non- conformal mapping, defined by a suitable 
function of p, 9 and T; 

jx = X(p, 0,T) 

, Y = 9 
T = T 


p= p(X,Y,T) 

0 = Y (5) 

T = T 


2 


will transform the region of interest in the right-hand side of the ''-plane 
bounded by 0 = b and p = c onto a rectangle, bounded by the lines: 

^ X = 0, corresponding to P = b (body) 

X = 1, corresponding to P = c (bow shock) 

j „ „ 

y = - corresponding to ® = - '2 (windward symmetry line) 

Y = corresponding to 9 = ^ (leeward symmetry line) 

An example of such a function, X(o, 9, T) will be discussed in Section XI. 
III. DERIVATIVES RELATED TO THE MAPPINGS 


Let 


dC _ iUJ 

8 = 31=0 = 


(7) 


be the complex derivative of C with respect to z (at t, t, T all constant); 
similarly, let 

CO = I = coi + 1CO3 (8) 

From (3) and (7) it follows that 

(9) 


Oi = ,1ZS_ = = C + iS 

U/g| 


where 


C = cos(9 - u)) , S = sin(9 - uj) (10) 

We introduce now the notations, ijf and f, for two analytic 
functions obtained by differentiating g and C with respect to t (that is, 
at constant x and y): 

1!/ = + ill's (U) 


f = i^=£, +if. 


( 12 ) 


Recalling that 


(14) 


St 


p ”5r 


+ i 


50 

W 


we obtain: 



1 P, = GC, 

p^=GS, 

Pt= Oh 



( 

0 = - 
X P 

0 = 5c, 

y P 

9^=fq 


(15) 


T = 0 , 

X ' 

7=0, 

y 

T 1 

t 



Conversely, noting that 






'‘t' - ®l>’ 

Vt = - (Vo 

1 >'8 

(16) 

we obtain; 






X = — C 
f ^p G » 

1 

p '< 

•'^0 = - G® 

, = - (0£i 



! 

i 

ly = i 3, 

1 ^P G » 

ve = 5=- 

Yt = - (® 



1 

1 


t0 = 0, 

t,= 1 



Between the two sets, (p 

, 0, t) and 

(X, Y,T), the following relations hold 


f 

P 

Py=-^p 

■ "t " ■ Xp 



< 

® ° ' 

9y= 

e^ = o 


(18) 


, = 0 , 
p- Px' 

T =0 

Py 

0 Px 

Tt= 1 

Pm 

, X = -^ 

T Px 



i 


Y0 = 1 , 

Y,= 0 


(19) 


T = 0 . 

p ’ 

T0 = 0 , 

T,= 1 




V 

By combining (17) and (18), we obtain: 
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^ ""x = Gir » 


s 

'^x ° o3r ’ 

*x = “ ■ 


’‘y“'5(°cX 7 + *) ■ x:"(”' 


>'y = b(-®^ + °) • 


ty = 0 


8 

" "G X 


‘t = ‘ 


I. (iS£,+ Cfe)g 


( 20 ) 


The following formulae are also obtained easily; 

G^. = GfC'i - 05 i fi + cca fa] 


r* - 


G0 = - Gcd 


9 > 


( 21 ) 


m^ = £®= , 


Wg = Gcoi , lU^ = t'a - CDx fa + Wa fi 


G^ " » Gy. = “(^<“''1 C-:Pn} , G^ = G 


(22) 


I U)^ = - ^ (^% - Ccpa ) , 


ujy = ^'(Ccpi + Scpa) , U)|, s i)(a 


IV. IMPORTANT UNIT VECTORS 
We begin this section by defining a p-line on a t=constant (physical) 
cross-sectional plane as a line along which 9=constant; similarly, a 
0-line will be a line on the t=ronstant plane along which psconstant. The 
unit vectors, i and j will be used to identify the tangents to a 6-line and 
to a p-line respectively. Note that 


I = C i - 5 j 

^ A A 

J = 3 i + C j 


i = C I + 3 J 

A A A 

j = . S I + C J 
By using (20 and (23) for any point, 
Q=xI+yJ+tK 

we obtain; 


(23) 


(24) 


5 




= " 5(p3r - §£3 j + K 

p 

The unit vector, N, normal to an X=constant surface, is impor- 
tant for the calculation of body and bow shock points. The body, indeed, 
is defined by X=0 and the bow shock by X=l. In general, 

1 


‘0 


N = Ni i + Na j + NaK = 


1Qy><Qt1 


Q^XQrj, 


where Q is a point on the surface. From (Z5) it follows that 
= ^ ' Na =-^Ni, N3 = Njd 


(26) 


(27) 


with 




■p ■ p 

In particular, at the body, from (18), 


X 


0 


oX, 


-Y 


X " “ *^T 
P 


(28) 


(Body) (29) 


Note also that we can write by or bQ, and b,p or b^, indifferently. 
Therefore, at the body (27) and (28) take on the form; 

, b„ 


Ni = V' 


Na = - -gi Ni , N 3 n Ni d 


(Body) (30) 


G^^T ■ Y 
Similarly, at the shock. 


b 

= - i(b,j,+ bYf 3 -bfi) , v = ,,/l + (-^)+ d^ (Body) (31) 
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I 


(32) 


^8 

o5T 




= - c. 


(Shock)< Ni = -• 


N, - - --iNi , Na :: Njd 


= -i(c„ + c„£,-c£i) 


c ^ 

V S /l+ +d^ 

V ' c / 


Let 

P(x,y,t) = 0 

define the geometry of the body in the physical space. The Image 
body in the (p, 0 ,t) space is 

P = b(8, T) 

To evaluate (30) and (31), that is the normal to the body, v/e need 
and b.,. At t= constant, 

‘ ^-*0 ‘=0 + + V’'o ‘’8 '' 0 ’ ‘ “ 

Con 55 ‘'»quently, and using (17)i 

bfl S.7 - C? 

9 X y 

"E" = 

X y 

Similarly, at 9=constant, 

? (x„b_ + x_) + 5^ (y*b + y } + ? = 0 

X p T T' y ' ' P T ' T t 

and 




b^ = - fa bg + bfl - G 
Therefore, at the body, 


_y . ^ 

X y 




t 


(Body) 

and (30, (31) can be replaced by the simpler expressions: 


^ " Cl? V 

X y 


(33) 

(34) 

(35) 
of the 

(36) 

bg/b 

(37) 

(38) 

(39) 

(40) 

(41) 
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1, N, = 21, Na = -± (Body) 

V = ^ 7 = +73+ 73 (Body) 

The above formulae are general. For any particular geometry, 
and ‘’'0 must be evaluated. 


( 42 ) 



V. EQUATIONS OF MOTION 

Having chosen a suitable reference length, 5^j.g£> the pressure, 
density and temperature of the free stream are chosen as reference 
pressure, density and temperature, respectively (Pyg£i Pj.g£> 

V.'-th p, p,7 measuring non-dimensional quantities, the equation of state 
is then 

p = py (43) 

The reference velocity, ^j.g£» is defined by 

"?e£=Pref/'’re£='*^re£ 

where R- is the gas constant divided by the molecular weight of air. The 
speed of sound in the free stream, in a non-dimensional form, is then 

= ./ Y (45) 

The logarithm of pressure is denoted by P; 

P = In p (46) 

A non-dimensional entropy, S (which is the difference between the local 
entropy and the free stream entropy divided by c^) is related to non- 
dimensional temperature and pressure by 

S = Y In y- (Y - 1)P, 7=exp(^P + is) (47) 

Euler's equations of motion in non-dimensional form are; 
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V . 7P + yV . V = 0 


i V(V ® )-VxVx^ + J7'7P=:0 


V . 7S = 0 


With k - K, let 


where 


and let 


V = w(x+ k) 


V = a i + Tij 


3 ^ 3 ^ 

= ' 55 ? ^ ^ 


Note that 


j V‘ ^P = w(x + k)* (’^iP + P^k) = wx‘ + wP^ 

\ V • V = • [w(x + k)] + vv^ = X* w + wV^ . X + 

j (52) 

) (V®) = w Vw (X® + 1) + iw3 Vi X® + i w= (^3)^ k 

I V xVx V = - w('!7i w X + w^)(X.+ k) + w^w(x® + 1) + w®[(X^ • x)k ~ X^ + ^ XJ 

therefore, (48) take the form: 

^ w(x* ^iP + Pp + y(X* w + wVj • X + = 0 

^ ^ w® X® + w » X + w^)X + X{. - w® X>c^i X X + P = 0 

w(Vj^ w *x + w^) + 7P^ = 0 

, b+S^.= 0 (53) 

The third of these equations can be used to simplify the first and second 
equation; finally, the following system is obtained; 
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( 64 ) 


(1 - — )P, + y. 7 iP + vVi.5=0 

{ v" - XX V, X y+-^(^iP - P^x) + X^ = 0 

w® 

s + s = 0 

\ ' 

The third of (53) is not needed; the above system is composed of four 
scalar equations for the two unknown scalar, P and S and the two- 
component unknown vector, X* Once P and S are determined, « is 
obtained from (47); the modulus of the velocity, q, is obtained from 

q” -7) (55) 

where 3"^ is the (non-dimensional) stagnation ‘cemperature, and w fol- 
lows from 

w^(l + a® + T'®) = q® (56) 

There are definite advantages in using (54) as a basic system of equa- 
tions (instead of (48) or of equations in divergence form). Not only it 
contains only four differential equations to be integrated, but it provides 
a clear separation of unknowns, S on one side and P and x on the 
other side, which is particularly welcome in problems where strong 
entropy gradients occur (Refs. 1,2). Another advantage of (54) stems 
from the fact that operates on the (x,y) plane only; therefore, it can 
be expressed in terms of p and 0 as independent variables, and using 

A A 

i and j as unit vectors. In particular, note that 
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7,P = o(Ppi“ + ipjj) 


V 


X.'?i P= G(aPp + ~Pg) 


5x7,. .(gn(.i-aj) 


( 57 ) 


7! _ ^ r^Pff 


^i’X=^[(-^) + (~) ] 


e 

I (X®) = G[(a Op + n Tip) i + I (a 7 q + r Tig) j ] 

X^ = [o^ + (uj^ - B^)v]i + [tTj. + (0^ - u)^)a] j 

These expressions can be substituted into (54)j in doing it, however, 
note that t also must be substituted by T, and that, for any function, f 

^t “ -p Pt ^ -6 ®t “ 

Using the notations; 


= I - — 
,2 


Bl. 


■A.1 = G + p fi , 


Bi = G o + p fj. , 


Ao - 5 4- f 

Aa - p H + h 


Bs = ^ -n + fa 


(59) 


B = ■p[Tl(l - Crtl) - CT Cpa] + fg - llfg 

and taking (21) and (22) into account, (54) become; 

Pt + ^1 Pp + Aa -“i)-''- Wa] = 0 

BsOgt — [- (jP^+ (O-o p£i)P - ofj Pg] - n D = 0 

71^+ Bi Tip+Ba rig+-^[-TiP^- pnfi Pp + (| - nfa)Pg] + ctD = 0 

V/ 

+ Bi Sp + Bg Sg = 0 
The final form of the equations of motion is obtained by expressing the 


(60) 
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derivatives in terms of X, Y, and T, considering that, for any function, 

X 

P X p 

= §T + 


and, consequently: 

$p + B:, = §,p + E + Bs 

+ Ai §p + Ag $g = $^ + C + As 

where 


(62) 


C = X,j. + As Xp + A^ Xg , E = X,^ + Bi Xp + Bs Xg 


With the additional notations; 

L = a(l - coi) + htpa 
F = - a X^ + (G - apfi) Xp - afg Xg 

H= - nX^ - pTlfi Xp + (|. nfg)Xg 


(64) 


the equations to be integrated at every grid point, except on the body and 
on the bow shock, are; 


/ + CP^ + AsPy+^[X a^ + i(rY + Xg^+ L) ] = 0 


a_ + Ea^+Bg 0Y+ — [- a P^, + F P^ - a fg Py] - tiD = 0 

w® 


(65) 


tIy+EHj^+Bs nY + ”[- TiPrr + HP^ + (^- nfg)pY] + aD = O 

w^ ^ 

\ S^ + ES^+Bg Sy = 0 
VI. CHARACTERISTIC EQUATION FOR BODY AND SHOCK POINTS 


Equations (65) are not suited for the numerical analysis of points 
on the body and on the bow shock. A characteristic equation in an (X, T) 
plane is needed; note that in the physical space the image of such a 
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plane is a surface almost parallel to the local velocity vector and almost 
normal to either the body or the bow shock. 

By multiplying the first three equations {65) by Mi, Us, and Wg, 
respectively, adding, and calling X the slope of the characteristic in 
the (X, T) plane (X = dX/dT), we obtain the compatibility equation: 

Fui - — (cMs + 'nU3)](P^+XPx) + Wa(cT^ + Xav) + U3{T^,p + Xr^) = 
w® 


= Uj + Ua Ra + W3 ^^3 

where 

*7 

( Rg = - Bg (Jy + CT fa Py + "nD 

w® 

. R3 = - Ba <p - (I - Tifa)PY ‘ 
w® ^ 

In turn, X is defined by 


(66) 


(67) 


C - X 

— (F + aX) 

— (H + -nX) 




w® 

w® 



YG„ 

E - X 

0 

= 0 

(68) 

yg„ 

iTp^e 

0 

E - X 




that is, 

X 

(E -X)(C -X) - ^ 2f_l (H+ nX) + X (F + aX)] = 0 
w® " ^ ^ 

or 

X 

X® - 2{X^+AiX +Aa Xg)X +E C-^ ^(FX +H-1) = 0 


which, after some manipulations, yields; 


X 


C + P 


aG 

WH 


(69) 
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( 70 ) 


/' ^9\ 


/ 3 % 


The lower sign and the upper sign must be used at body points 
and bow shock points, respectively. 

From (68), we obtain 


Ui = E - X 

11 — 7G yr 

11 - 7G v 

- “ iTp-^e 

Wi - — (CJUs + t1|i3) = + P ^ 
w® 


The compatibility equation (66) is then 


+ (P,^ + XP^) = +-^('r'T+^ V 


X, 


(71) 


(72) 


+ [-^ -T-- 

r K 


-i?Rj 


(73) 


VII. EQUATIONS FOR BODY POINTS 


At body points, the boundary condition, 

V • N = 0 (74) 

yields 

b 

cr-'p-r- +d=0 (75) 

D 

if (49), (50), (30) and (31) are taken into account. On the other hand, (63,' 
and (29) give 

E = X^(-b^ + Bi - Bs b^) (76) 

and, using (59) and (31), it is easy to see that 

E = 0 (77) 

It is also important to note that (75) is identically satisfied for any 
value of Y and T; in particular, 
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(78) 


^ Y 

a„ - - 5 “ = 


The compatibility equation (73) becomes 

b 

jf. 

T 

with 

aG 


(79) 


X :r C - P 


WH 


(80) 


Y 

It remains to provide explicit expressions for (-^) and d^. 


The former is easy to write: 

>b-,\ , / b. 


(-"I = 

To get the latter, we may start noting that, by (20), 


( 81 ) 


+ [sb^.b(*£,+C£a)]dy j 


(82) 


Using d in the form (41), we obtain 

■>T = cTT^l^’tt-’^'Vx-^ + V ^ 

X y ( ^ 

+ i<rf?b_-b(Cfi -§fp)T? . -d(f 5^ +3 ? + C7 +3? )1 + 

Gi- T ‘’ji.xt :xx xy xx xy'J 

+ [Sb^ - b(®fi -d(fy?^+ + ssr^) 

= -n) -d(C3'^, + S7 ) + 

X yl ^ ^ 

+ 1 < [c bT - b (C fi - 3fa )][[y^ - d(S - C )^[S ( 1 -toi ) + C d (C 3 ?^^) ]] + 

+ [S b,^- b(3 £, + n fs )]IT!?y^+ d(3 7^- G ?y)§[c(l ) -3 r 03 ]-d(G S^^y)]]) 

That is, 
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* ■'yt> 4^' '’4- '=<'’ L%r ’xx-' ■' 'xyl] ^ 


t>Y'“ 


+ L® - b('? + f’ f. ) ; d(f’ 3 9^^)j ) j + 

b„ 


+ d -gi,- •'£. - ^-5 + fi ra + fg ri ^ (83) 

Having used (79) to determine P^, the last of (65) yields S and 
from updated values of P and S, updated values of S' and q® are obtained 
from (47) and (55), 

Now, let 

b 

V s w (r + a) (S4) 


The third of (53) in the (X, Y, T) frame, with E=0, reads: 

w^ 4 Bg Wy +|[Py + (X^+oXpfi + Xgfg)P^ + faPY] = 0 

(85) 

If (85) is multiplied by and added to the third of (65) and the second of 

w® 


• 1 

(65) mviltiplied by the following equation is obtained 


r,by 


+ +B, (^) ]aw + ^Py + Dw(o-^r) = 


SG 


( 86 ) 


which can be used to update v. Note that v is the velocity component 
tangent to the body in the cross-sectional plane. Its Lagrangean derivative 
expressed in the (X, Y, T) frame by v^ + Bg Vy» depends on the geometry 
of t-*e body and on the Y-derivative of P only. After updating v, updated 
values of w, a and v at the body can be obtained by solving (84), (56) and 
(75) . One obtains; 


w = 



with V defined by (31), 


(87) 
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(v/w){b /b) - d 

a = 1 (88) 

1 + (bY/b)® 

" b^ 

m 

VIII. BOW SHOCK CALCULATION 

Let 

7i^ 71^ i + % ] + Tla k (90) 

be the unit vector normal to the bow shock surface. The values of 

are the same as the values of Nj , Na , No defined by (33). The 
velocity component normal to the shock in front of it, u^, is 

^co = ^ ^3 (91) 

where 

u = V S sinCX 

CO 00 

V = V C sinoc 

00 CO 

W = V cos a (92) 

00 CO 

If we denote by u the corresponding velocity component behind the shock, 
the velocity vector, V behind the shock is; 

V = V^+(u-uJ^ (93) 

CO CD 

The Rankine-Hugoniot conditions provide the increment in P and the ratio 
u/u^ across the shock (here P is the logarithm of pressure behind the 
shock; let us keep in mind that P^^ = 0): 


p = In V J 1 + 2 ‘) 

(94) 

v-1^ . 2 y 1 

" = V + 1 + V + 1 - 

(95) 


Since (94) and (95) are identically satisfied at any T, 
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(96) 



2 u. 


" 2 


TTI '^toT 


u. 


• ’ill 

v+i 


2v 

v+1 



(97) 


In turn, from (92), (9)> (18) and (21): 




V „ s u OJm 
obT « T 


(98) 


'^«T = ° 


and 


Wm = ilfa - f» "-'^1 + + Go5a 


(99) 


Therefore, 

"mT “ “ **■ ^ ^ (100) 

In all the above expressions, we will now separate from the 

rest, letting first d^ in the form: 

d^ = Cl + Cs (101) 

The values of Ci and Cg will be computed later on. It follows that 

V(j) “ G3 + C4 (102) 

C3=i[^(^) +dCil, C4=^C2 (103) 

V C C fji -J V 

From (33), 



^6 '^TT ’ 

Cs 

- ^ c 

^3 ) 

Cs = - — C 4 

( ^3 T “ 

Cg t 

C 7 


Y 

Ce = - Cs 

7^3 T = + 

\ ^ 

Gl 0 C pp t 

Cg 

“ Cg d + ^1 Cl , 

Cl 0 — Cg d + 71^ Cg 


( 104 ) 
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Then 


i ^TT 

(Cn = + v^ Cg + C? + w^ Cq 

>Ci3 s Co + Cg + Cio 

From (96) and (97) 


(105) 

(106) 


2vico 

Fiji ® ^13 4 ^'P'P * ^13 “ V - 1 

- ”2~ 


2u 

u^V 


u^ = C16 + Cia , CiB = ~)cii , Cio= (^-:^~)Ci 8 


Then, 


(107) 


(u. - “ ^1’!' ***^18 ^tT ' Ci7 s CiB - Cii , Cje = Cjo - Cja 

(108) 


With 


^ A ^ 

V = ui + vj + wk 


(109) 


it follows from (93) that 


Therefore, 



Cg 0 ^'Ji'Ji 
■f Cga ^ipip 

+ Cp4 


and 


u = + (u-u„)?^i 

V = -f (u-uj?la (110) 

w = w^ + (u-u ^)?^3 


Cio = - U)rn+ Ci7??i + (u-u )Cb , C = Cia??i + (u-u )Ce 

X J3 Q W 

Cai = + Ci7?Zo +(u-Uj,j)C7 , Cop = Cl e?^p + (u-Ujjj)Ce 

Cq 3 = Ci7 ?^3 + (u-U^)Cg , Cp ^ "CigJ^a + (u-U^)Ck 


(111) 
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= Cob + Cc 3 » Cac = ;^(Cio - t Con) , Coo = ^ (Cqo * 7 Cy^) 

~ ^TT ' ^ C03) , Cyo = “ (Coo ■ C34) 

(112) 


By substitution into the compatibility equation (73), we obtain 

'''^0 r *^v 

" (C^ 3 + Cj 4 C^,j, 4 * \Pj^) 2 jr ' I^Cg g 4 *Coo + ^ (Cji 7 + Cg q ^PP 

Tn® ^ V “* C y ^ 

_,^(cr-in--i) + a0 Ri -R, f^Rgr 

1 * W C w • 


which allows 

Ctt 

D 7 


to be computed, as follows; 

~ D 7 (Dg + De + 7T) 

YwX“ cT^ 

Ci4 +— p-^(Caa - ~Cge) 


(113) 

(114) 




VWX. 

c„ 


Dg = 

Ci3 + 

-^-^(Cas - 

— Cg7) 

(115) 

Da = 

VwXp 

fc)a 

{-H. 4 ^ 3 } 

r* sX ^ C.. “1 

- 0^1 

(116) 



YwX„ 

*1 


7T = 

’'[px 

Ta ^^X ' 

- -r ’■xO 

(117) 


Having computed ^«p ^.nd c are obtained by successive 

integrations. Once a new shock geometry, p=c(9, t), is obtained, the 
A ~ 

new ?l and are evaluated; then, (94) and (95) can be applied to compute 
P and u. The three velocity components follow from (93); therefore 
a(=u/w) and r(=v/w) are made known, Finally, S is given by 

S = P - V In— (118) 

u 

To complete the information needed for coding the shock calculation, 
we need to evaluate C^ and Cg . From (34), 
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d rji ^ ^ Grji Cijnip + c 

S -d['’'i +i(C"'i +S«)3)^Cc,^-c(Cfi - 3fa)| +i(3'45i - Ccpa )(^ - c(S +(’fa)}' 

'* G *^TT ' G YT ) " G Y T^ 

■ ^ T *1 1 1 

= - d|^ C'l +-^C0i " 'Sifi + JCs^a j " ^‘''j'T ” ^ (cyiji^s ■ " 

- ■^[c^ Imag (f^) - c Real (f^)j (119) 


Now, 

£t- f^. + fjjXT+ y^ 


ot 

= - c(C f, . 8 )]+l -^(f)[s = J . c(? £, . C £. ) ; 


.S!^c,^l|„.£,T2i,eT.cf): 

Therefore, 

Cl = > dfilfi +-^cpi - tOi £i + CDs fa^ -■^ (cy^fa - c^ fi ) - 
- Cy Imag -^Imag'[(\i! - f)(Cy - cf); - 


Real -^^i2£i-Real((>l' -£)(c„-cf)^' 

Bt= >. T 




Ca = -^ 


(120) 


( 121 ) 


IX. GENERAL OUTLINE OF ONE INTEGRATION STEP 


The equations obtained in the preceding sections are used to 
proceed from a station, t, to a station, t + At, using a predictor-corrector 
integration scheme, as follows. 
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Predictor stage 

Given, original valueo of P, C, S, q® , JT, w, c, c,^,, b, b^, 

7'» g» f, and it< 

M all value of 6, compute Cy, Cy^, byy, byy. 

At all value of 9 and 0 

Compute X- and Y- derivatives of P,0, r, S (see next section) 
Compute p = mod C, X^, Xq, X^. (see Section XI) 

Compute C, H, Aj.Ab, , Bg , D, E, C, L, F, H 

For all points except body and shock points, determine Py, Oy, 

riy and Sy from (65) 

For body points, compute , Rg , Rg from (67), \ from (80) and 
dy from (83). Then, xise (79) to obtain Py, (86) to obtain Vy, and 
the last of (65) to obtain Sy. 

For shock points, compute Ri,Rg,R3 from (67), evaluate 7Zi , Tig 
and Tig and compute all coefficients (C^ through Cga). 

Compute D7,Da,Eb and FI, and determine Cyy from (113). 

Update P, 0, n and S at all interior points, P, v and S at all body 
points, as well as Cy and c, using the following rvile (where ® is an 
arbitrarv function): 

$(t At) = $(t) + $y At (122) 

Move to the station defined by t + At; compute the geometry of the 
body; determine the basic parameters for a new mapping. Determine new 
values of b,by,by. From the updated values of c, determine new values 
of Cy; using the updated values of Cy, evaluate the new values of TIi,TIb, 
Tig, the corresponding u^^ from (91) and updated values of P and u behind 
the shock from (94) and (95); then use (110) to get updated value s of u, v, 
and w and update a and ri accordingly. Use (118) to get an updated value 
of S behind the shock. At the body, use (87), (38) and (89) with the 
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updated values of , v and b, At all points, updated values of S', q® 
and w must be computed. 

Symmetry conditions are now imposed where necessary. Finally, 
the original values of P, a, v S, v, c and are temporarily saved. 

A new computational grid (including values of g, f, cp and it’) is generated. 

Corrector stage 

The computation is restarted as at the beginning of the predictor 
stage. The geometry and the grid, however, are now those of station 
t + At; and all the variables have their predicted value at t + At as well. 
The updating in the corrector stage is performed using the following rule 
in lieu of (122): 

5(t + At) = ■^'l_$(t) + T(t + At) f $^At] (123) 

values $(t + At) is the left-hand side of (122) and is the T-derivative 
computed In the corrector stage. A new evaluation of the body geometry 
is not necessary; the computational grid, however, has to be re-evaluated 
since the bow shock location may have changed slightly. The values at 
the shock themselves have to be recomputed as at the end of the predictor 
stage. Symmetry conditions are imposed again. Finally, the updated 
values of P,cr, ri,S,v, c and c^ are stored as initial values for a new step, 

, t 

X. DISCRETIZATION OF X- AND Y- DERIVATIVES 

(3) 

In principle, the MacCormack scheme' ' for integrating the 
equations of motion at all interior point is adopted. This is reflected by the 
use of (122) and (123) at the predictor and corrector level, respectively. 

In addition, whenever possible, the X- and Y- derivatives should be 
approximated by 2-point differences taken, for example, forwards at the 
predictor level and backwards at the corrector level. 

There are several exceptions to the above rule. 
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1) At body boundary points, where X=0, backwards approxima- 


tions cannot be taken and are replaced by the approximation: 

-E$o + - $3 


$ 


X 


AX 


(124) 


which maintains a second order accuracy if the equations are linear 

(4) 

throughout the integration step' \ Here, §. = $(iAX, Y) 

2) Similarly, at bow shock boundary point*’, where X=l, forward 
approximations cannot be taken and are replaced by the approximation: 

2$o - 3§i + $3 


6 PS 

■X 


AX 


(125) 


where = §(1 - i AX, Y) . 

3) The last of (65) simply expresses the fact that entropy is 
transported unaltered by the moving particles. That means that no 
entropy signals travel backwards along streamlines. Consequently, no 
approximations to X- or Y- derivatives are allowed which imply forward 
differences of S along a streamline. 

In most of the cases the rule does not have to be strictly enforced; 
in other words, the violation of the physical principle produced by alter- 
nating backward and forward differences (as in the MacCormack scheme) 
does not produce sizable errors, so long as the entropy distribution is 
smooth and essentially flat. Difficulties and inaccuracies appear, however, 
when entropy layers tend to build-up^^^. It is therefore advisable to adopt 
the following scheme systematically: 


So 0 - Si 0 

iPredictor level: S^sssgnc — 


2Soo- SSio + Sso 


Correction level: S^sisgnO' S'X" 


So o “ So 1 
S^^sgari ^ — 

2Sqo ~ 3So 1 + Sq a 

Sy^^sguT 


( 126 ) 
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where 


Soo=S(X,Y), Sio - S(X- iXsgnC, Y) 

S,o = S(X - ZilXsgnC, Y) (127) 

So 1 = S(X, Y - AY sgn r) , S 09 - S(X, Y - 2 AY sgn v) 

At body and bow shock boundary points, the rules given at (1) and (2) above 
should be used for approximating as well as the X-derivatives of the 
other physical parameters, 

4) Similarly, the terms v^ + Bg in (86) express the Lagrangean 
derivative of v. Forward differencing of v is thus forbidden and Vy 
should be approximated as follows: 


Predictor level: 


Corrector level; 


^ Vo - Vi 

Vy wsgnv — ^ 

2vg - 3vi + Vs 
Vy«sgnv 2 Ty 


(128) 


with 

Vo = v(0, Y) 

Vj^ = v(0, Y - AY sgn v) (129) 

Vg = v(0, Y - 2 AY sgnv) 

5) Y-dorivatives at the bow shock boundary are conveniently 
approximated by centered differences . 


XI. COORDINATE NORMALIZATION AND GRID STRETCHING 

ALONG p-LINES 


The object of the transformation (5) is twofold. It defines a 
variable X which is constant (equal to zero) along the body and also con- 
stant (equal to 1) along the bow shock. In addition, it provides a stretching 
of coordinates according to which evenly spaced grid points on the X-axis 
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correspond to unevenly spaced points on the p-lines. The latter property 
is used to accumulate 6 -lines in the vicinity of the body where a stronger 
resolution is needed. 


The values of the derivatives, Xp, Xq, and X.^. depend on the 
choice of the stretching function X(p, 8,t). in this section we give an 
example of such a function. If the definition of X(p, 0, T) is changed, the 
definitions of Xp, Xg and X.^ must be changed accordingly. The rest of 
the program does not need to be altered. 


Let 


p = c + tanh [a(X- 1)] 


(130) 


with 


c-b 

tanhO. 


(131) 


Obviously, X as defined by (130) satisfied conditions (6). Different 
values of cx provide different degrees of accumulation of 6-lines near 
the body. To give an idea of the effect of Ot, let b=0, c=l. Fig. 1 plots 
X vs. 0 for various values of Ct. Clearly, strong stretching effects begin 
to appear for values of Q, larger than 2. 

It is easily proven that; 

1 


X = — 

- (p-c)®] 


X, 


“ " {tanha ^ “^yI 


X. 


I tanha "''Y 
r p-c 


" 1 tanha ^^T " ^T^ ^ '^tI 


(132) 

(133) 

(134) 
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XII, EXPLICIT COMPUTATION OF TERMS RELATED TO THE MAPPINGS 
The mapping of the z-plane onto the ^-plane, mentioned in 
Section II, is performed according to the general scheme exposed in Ref. 

5. In the z-plane, the "hinge -points" are denoted by h^ ^ (-1=1 through J). 

In addition, h ^ ^.and j ^ire the affixes of the lower and upper 

intersection of the cross-sectional contour with the y-axis. 

Let 


= z 


(135) 


and a sequence of J mappings be used, each defined by the equation; 


■j+l + ^ \: 


6 . 

z. - h,. \ '^ 

} 11 


1 

\z. + h.. 
\ J JJ 


(j=l,2,.. ., J) 


(136) 


or its inverse: 


z . - h.. /z, , - 1 

11 - / . 1+1 ^ ■' 

z. + h*. r +1 + * 

J JJ 


(137) 


The mappings are automatically performed in order of increasing 6,. 


With 


^ 2 ^^J+1, J+l + ^J+2,J+1^ 


the £ -plane is defined as 


^ " ^J+1 " ^ 


(138) 


(139) 


To obtain the derivatives used in the computation of the flow field, 
the following definitions are used: 


h. . = Ct. + i P . 
JJ J J 


(140) 
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h.. - z. 

k, = -U i- 

a. 

J 

ko = - 1 

k 3 = - a. k, = Z. - h.. 

k 4 = z. + h.. 

J JJ 

ks = log (ks/k^) 

k, = k. 

k7 = z. - h.. 

J JJ 

ka = z . - i P . = ka + a , 

J J J 

kg = z. - i p. = k? + a. 

J J J 

4 4 

ki 0 = 6./6. +ct,/c(,, 

J J J J 

kii= l/{kak4) 


(141) 


where dots mean partial differentiations with respect to t at constant 
X and y, that is, at constant z). 

We also define 

d z.j^^ 

(142) 


d z * , - 

g,- = I's ^11 


’ j d z . j J 


and we note that 


h„ + — = gi(ki ke - h. ) 


® j+ii 4. 'Jim-* 


Therefore, 


M Sh?', M “) 

JJ 


= g^.(k7 + ki ko) + 2 ks ka 6^ 


JJ 


(143) 


( 144 ) 


and 


4 4 


1 


- ' “^J+l ■ 2 ^^J+1, J+1 ^J+2, J+1^ = ^'j+1 “ ^ 
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The t- derivatives of the h^^ can be obtained as particular cases of (144), 

with h^ h^j in lieu of respectively. In particular, if -tsj, 

« 

h. . , , = 0, as evident from (136) and h. s 0. 

J/J+1 

Second derivatives are also needed to evaluate 5° log {!/5t^ which 
is used in (121). From (142), it follows that 


ka ^j+1 ^j+1 ■ ^ ftj)J 

and, from (144), 


(146) 


- h. • - — (k7 + kj ka ) + kj ] + 


= g.(k7+ka ke) 

L^j+1 ^j + 1 1 j ke 6j (147) 

As above, the second t- derivatives of the h^ can be obtained as particular 

i 

cases of (147), with h^ , h^j in lieu of z^ respectively. 

At every new cross-section, the geometry defines the values of 

• • M 

h^l, h^j and h^j^ (-6=1 through J+2) as well as of 6^, ty and 6 j (j=l 

through J). According to the procedure explained in Ref. 5, the mappings 

are performed in order of increasing 6.. First, the values of h, 

J •('.J+l 

(j=l through J) and found by repeated applications of the mapping routine, 
which also provides the values of ^h^ 

Then, E is obtained from (138). Repeated applications of (144) 

• •• 

and (147) allow all the values of h,. and h , . to be determined. It is 

-t-j 

• .. 

possible, thus, to evaluate E and E as well. 

Grid points on the contour of the body are determined next. For 

i0 

each value of Y (that is, of 6), a guess of b is made and C = b is 
used to get through (139); then, the sequence of mappings (137) with 
decreasing values of j Is applied to obtain a value z which should match 
a body point in the physical cross-section. Matching within a prescribed 
tolerance is obtained by a trial-and-error procedure. 
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Knowing c{6) and having just determined b‘(9), values of o are 
obtained at every grid point outside the body by applying (130) or any 
equivalent expression if another stretching function Is used. For such 
points the mappings are applied again as for the body points but no trial- 
and-error procedure is necessary. The entire grid is thus obtained both 
in the ^-plane and in the z-plane (the physical cross-sectional plane). In 
addition, at this stage all values of g. have been determined by the map- 

J 

ping routine. Consequently, it is possible to apply (144), in order of 
increasing j, to obtain all the values of and, similarly, (146) and (147) 

M 

to obtain all the values of The latter calculations are actually per- 

formed only at shock points, the sole points where second t-derivatives 
are needed. The functions g,cp, f and '■1' necessary in the aerodynamical 
part of the code are now determined at each grid point as follows: 



J 

- 2i V 

L. 


g 


J 


r 


j=i “j+i 


1 


^j+i 




f = 


1 

C 'Sf ■ C , ^J+l 



J 



where the right hand side is obtained from (246). 


(g,= l) 


(149) 


(150) 


( 151 ) 
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Finally, the term ^’^log'/St^ which appears in (121), is com- 
puted as follows: 


loR ^ 


c£ 1 


= = J+1 ■ 



(152) 
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